/**
 * @author ahammer
 * @date   2010 01
 * @update 2010 01 added progress bar support
 * @update 2011 01 added 16bit option
 */

/*
 * The implemented method uses IFR Inverse Fourier Reconstruction  with 1/(f^(b/2)), filtering and b = (8-2D)
 * see book: Turner, Blackledge, Andrews, Fractal Geometry in Digital imaging, S.172
 * see also Anguiano etal. J. Microscopy 1993, Vol 172, 223- 232 (FFT Dimension etc.)
 * and book Voss R.F. 1985
 */

/*
 * This file is part of Iqm.
 * Copyright (c) 2010-2011 Helmut Ahammer
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package op;

import java.awt.Point;
import java.awt.image.ColorModel;
import java.awt.image.DataBuffer;
import java.awt.image.DataBufferFloat;
import java.awt.image.SampleModel;
import java.awt.image.WritableRaster;
import java.awt.image.renderable.ParameterBlock;
import java.util.Random;
import javax.media.jai.JAI;
import javax.media.jai.ParameterBlockJAI;
import javax.media.jai.PlanarImage;
import javax.media.jai.RasterFactory;
import javax.media.jai.RenderedOp;
import javax.media.jai.TiledImage;
import javax.media.jai.operator.IDFTDescriptor;

import main.Manager;


/**
 * This is the main image processing class
 * There exist two approaches:
 * A user defined JAI operator is just called 
 * or:
 * The actual processing is implemented in this class 
 */
public class IqmCreateFracSurfOperator{ 



public IqmCreateFracSurfOperator() {
	  //WARNING: Don't declare fields here
	  //Fields declared here aren't thread safe!   
  }
  /**
   * 
   * @param ParametrBlockJAI pb
   */
  public Object run(ParameterBlockJAI pbJAI){
	    
	    Object ob = null;
	    //ob = JAI.create("IqmCrop".toLowerCase(), pb, null); 
	    
	    //PlanarImage pi = (PlanarImage) pbJAI.getSource(0);
	  	int   width  = pbJAI.getIntParameter("Width");
	  	int   height = pbJAI.getIntParameter("Height");
	  	float fracD  = pbJAI.getFloatParameter("FracD");	
	  	int   outbit = pbJAI.getIntParameter("OutBit");

	  	int greyValueMax = 0;
	  	if (outbit == 0){ //8bit
	  		greyValueMax = 255;
	  	}
	  	if (outbit == 1){ //16bit
	  		greyValueMax = 65535;
	  	}

	  	//DFT image width and height is only possible in 2^n steps.
	  	int dftWidth = 1;
	  	int dftHeight = 1;  	
	  	while (dftWidth  < width)  dftWidth  *= 2;
	  	while (dftHeight < height) dftHeight *= 2;
	  	if (dftWidth > dftHeight) dftHeight = dftWidth;
	  	if (dftWidth < dftHeight) dftWidth  = dftHeight;
	    
		//create real and imganiray part
	    SampleModel smN = RasterFactory.createBandedSampleModel(DataBuffer.TYPE_FLOAT, dftWidth, dftHeight, 1);  //1 banded sample model
	    SampleModel smM = RasterFactory.createBandedSampleModel(DataBuffer.TYPE_FLOAT, dftWidth, dftHeight, 1);  //1 banded sample model
	    float[] floatArrN = new float[dftWidth*dftHeight];
	    float[] floatArrM = new float[dftWidth*dftHeight];
	    DataBufferFloat dbN = new DataBufferFloat(floatArrN, dftWidth*dftHeight);    
	    DataBufferFloat dbM = new DataBufferFloat(floatArrM, dftWidth*dftHeight);    
	    WritableRaster wrN = RasterFactory.createWritableRaster(smN, dbN, new Point(0,0));
	    WritableRaster wrM = RasterFactory.createWritableRaster(smM, dbM, new Point(0,0));
    
	    //generate random pixel values
	    Random generator = new Random();		
	    //Random generator2 = new Random();		
	    float b = 8.0f -(2.0f * fracD);
		for ( int k1 = 0; k1 < dftWidth; k1++) {
			int proz = (k1+1)*100 / dftWidth;
   			Manager.setProgressBar(proz);
			for ( int k2 = 0; k2 < dftHeight; k2++) {
				float k = (float) Math.sqrt( Math.pow(Math.abs(k1-dftWidth/2) + 1, 2) + Math.pow(Math.abs(k2-dftHeight/2) + 1, 2) );
				//k[1,.......]
//				if (k == 0.0){ 
//					k = Float.MIN_VALUE;
//					System.out.println("k = Float.MIN_VALUE: " + k);
//				}
				//k = k/dftWidth;  //not necessary
				double g = generator.nextGaussian();
				double u = generator.nextFloat();
				double n = g*Math.cos(2*Math.PI*u);
				double m = g*Math.sin(2*Math.PI*u);
				n = n*Math.pow(k, -b/2);
				m = m*Math.pow(k, -b/2);
				wrN.setSample(k1, k2, 0, (float)n); 
				wrM.setSample(k1, k2, 0, (float)m); 		
			}
		}
	
		//create tiled images
	    ColorModel cmN = PlanarImage.createColorModel(smN); //compatible color model
	    ColorModel cmM = PlanarImage.createColorModel(smM); //compatible color model
	    TiledImage tiN = new TiledImage(0, 0, dftWidth, dftHeight, 0, 0, smN, cmN);
	    TiledImage tiM = new TiledImage(0, 0, dftWidth, dftHeight, 0, 0, smM, cmM);
	    tiN.setData(wrN);
	    tiM.setData(wrM);

	    ParameterBlock pb = new ParameterBlock();			
	    //Periodic shift: center is shifted to the corners of the image---------- 
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(tiN);
		//pb.add(shiftX); //default: Width/2
		//pb.add(shiftY);	//default: Height/2
		PlanarImage piN = JAI.create("PeriodicShift", pb);
		
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(tiM);
		//pb.add(shiftX); //default: Width/2
		//pb.add(shiftY);	//default: Height/2
		PlanarImage piM = JAI.create("PeriodicShift", pb);
	    
	    //create IDFT image

		//Reconstruct--------------------------------------------------------------
		pb.removeSources();
		pb.removeParameters();	
		pb.setSource(piN, 0);
		pb.setSource(piM, 1);
	    PlanarImage piDftNew = JAI.create("bandmerge", pb, null);		
	  
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(piDftNew);	
		pb.add(IDFTDescriptor.SCALING_UNITARY);	
		pb.add(IDFTDescriptor.COMPLEX_TO_REAL);
		ob = JAI.create("IDFT", pb);
			
		//Normalize grey values---------------------------------------------------
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(ob);
		RenderedOp extrema = JAI.create("extrema", pb);
		double[] allMin = (double[])extrema.getProperty("minimum");
		double[] allMax = (double[])extrema.getProperty("maximum");
		double[] scale  = new double[allMin.length];
		double[] offset = new double[allMin.length];
		for (int i = 0; i < allMin.length; i++){
			scale[i]  = (float)greyValueMax/(allMax[i]-allMin[i]);
			offset[i] = ((float)greyValueMax * allMin[i])/(allMin[i]-allMax[i]);
			//System.out.println("IqmCreateFracSurfOperator: allMin[i] allMax[i] " + String.valueOf(i)+ ":  " + allMin[i] + "  "+allMax[i] );
		}
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(ob);
		pb.add(scale);
		pb.add(offset);
		ob = JAI.create("rescale", pb);
		 
		//Change float to byte----------------------------------------------------
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(ob);
		if (outbit == 0) pb.add(DataBuffer.TYPE_BYTE);
		if (outbit == 1) pb.add(DataBuffer.TYPE_USHORT);
		ob = JAI.create("format", pb);
		  
		//Crop to original size----------------------------------------------------
		pb.removeSources();
		pb.removeParameters();
		pb.addSource(ob);
		pb.add(0f);
		pb.add(0f);
		pb.add( (float)width);
		pb.add( (float)height);
		ob = JAI.create("Crop", pb);
	    //--------------------------------------------------------------------------	
		
		//ob = piN;	  	  	
	  	//System.out.println("IqmCreateFracSurf ob.getClass()" + ob.getClass());
		((PlanarImage) ob).setProperty("image_name", "Fractal D="+String.valueOf(fracD));
		((PlanarImage) ob).setProperty("file_name", " ");
			
  return ob;
  }
}
